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ABSTRACT 

Finite  amplitude  solitary  waves  in  water  of  arbitrary  uniform  depth  are 
considered.  A  numerical  scheme  based  on  collocation  is  presented  to  calculate 
the  highest  solitary  wave.  It  is  found  that  the  ratio  of  the  amplitude  of  the 
wave  versus  the  depth  is  0.83322.  This  value  is  about  0.006  higher  than  the 
values  obtained  by  most  previous  investigators.  In  addition  another  numerical 
scheme  based  on  an  integro-dif ferential  formulation  is  derived  to  compute 


solitary  waves  of  arbitrary  amplitude.  These  calculations  show  that  the 

\ 

results  of  Longuet-Higgins  and  Fenton  (1974)  are  not  accurate  for  very  steep 


waves.  Graphs  and  tables  of  the  results  ar«  included. 
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SIGNIFICANCE  AND  EXPLANATION 


Since  the  tine  of  Scott  Russel  (1845)  many  approximate  solutions  for 
solitary  waves  have  been  obtained.  Most  of  these  calculations  are  in  good 
agreement  for  relatively  small  values  of  the  wave  height.  However,  some 
discrepancies  between  these  calculations  appear  as  the  wave  of  maximum  height 
is  approached.  In  the  present  paper  we  present  accurate  numerical  methods  to 
compute  steep  solitary  waves.  We  show  that  the  ratio  of  the  amplitude  of  the 
highest  solitary  wave  versus  the  depth  is  0.83322.  This  value  is  about  0.006 
higher  than  the  values  obtained  by  most  previous  investigators. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  of  this  report. 
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*  ** 

J.  K.  Hunter  and  J. -M.  Vanden-Broeck 

1.  Introduction 

Since  the  time  of  Scott  Russel  (1845)  many  approximate  solutions  for 
solitary  waves  have  been  obtained.  Solutions  in  the  form  of  an  expansion  in 
powers  of  the  wave  amplitude  were  derived  by  Rayleigh  (1876),  Korteweg  and 
de  Vries  (1895),  Keller  (1948),  Laitone  (1960),  Fenton  (1972),  Longuet-Higgins 
and  Fenton  (1974),  Witting  (1975)  and  others.  On  the  other  hand  direct 
numerical  calculations  were  attempted  by  Yamada  (1957),  Lenau  (1966),  Yamada, 
et  al  ( 1 968 ) ,  Byatt-Smith  (1971)  and  Byatt-Smith  and  Longuet-Higgins  ( 1 976 ) . 

A  review  of  these  investigations  can  be  found  in  Miles  (1980). 

Most  of  these  calculations  are  in  good  agreement  for  relatively  small 
values  of  the  wave  height 


Here  A  is  the  elevation  of  the  crest  of  the  wave  measured  from  the 

undisturbed  level  of  the  free  surface  and  H  is  the  undisturbed  depth. 

However,  some  discrepancies  appear  as  the  wave  of  maximum  height  is 

approached.  For  example  the  following  numerical  values  for  the  maximum 

amplitude  a  have  been  obtained!  0.827  ±  0.008  (Yamada  (1957)),  0.827 
max 

(Lenau  (1966)),  0.8262  (Yamada,  et  al.  (1968)),  0.827  (Longuet-Higgins  and 
Fenton  (1974)),  0.8332  (witting  and  Bergin,  unpublished  work  mentioned  by 
Witting  (1975))  and  0.8332  (Fox,  unpublished  dissertation  mentioned  by 
Schwarts  and  Fenton  ( 1 982 ) ) . 
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The  Boat  reliable  calculations  for  steep  solitary  waves  are  those  of 
Longuet  Higgins  and  Fenton  (1974)  and  Byatt-Smith  and  Longuet-Higgins 
(1976).  Both  calculations  predict  that  the  highest  solitary  wave  is  not  the 
fastest.  However  the  results  predicted  by  these  calculations  do  not  agree  for 
very  steep  waves  (see  Figure  3).  It  is  worthwhile  mentioning  that  Witting 
(1975)  has  suggested  that  the  method  used  by  Longue t-Higgins  and  Fenton  (1974) 
is  defective  because  the  assumed  expansion  is  incomplete. 

In  this  paper  we  present  a  numerical  scheme  based  on  collocation  to 

compute  the  solitary  wave  of  maximum  height.  The  method  is  akin  to  that  of 

Lenau  (1966).  However  our  results  are  more  accurate  since  we  retain  up  to  100 

terms  in  the  power  expansion  whereas  Lenau  retained  only  9  terms.  It  is  found 

that  a  =*  0.83322.  This  value  is  about  0.006  higher  than  the  values 
n unx 

obtained  by  Yamada  (1957),  Lenau  (1966),  Yamada,  et  al  (1968)  and  Longuet- 
Higgins  and  Fenton  (1974).  On  the  other  hand  it  agrees  with  the  values 
mentioned  by  Witting  (1975)  and  Schwartz  and  Fenton  (1982).  We  also  show  that 
Yamada' s  (1957)  scheme  yields  the  value  0.833  when  a  sufficiently  large 
number  of  mesh  points  is  used. 

In  addition  we  present  another  numerical  scheme  based  on  an  integro- 
differential  equation  to  compute  solitary  waves  of  arbitrary  amplitude.  The 
method  is  similar  in  philosophy  if  not  in  details  to  the  scheme  derived  by 
Vanden-Broeck  and  Schwartz  (1979). 

Following  Longuet-Higgins  and  Fenton  (1974)  we  introduce  the  parameter 


2 


Here  qc  is  the  velocity  at  the  crest  of  the  wave  and  g  is  the  acceleration 
of  gravity.  The  parameter  u  varies  between  0  and  1  as  the  wave 
amplitude  varies  from  zero  to  its  maximum  value. 


The  numerical  solutions  of  our  intsgro-differsntial  equation  differ  from 
the  results  of  Longuet-Higgins  and  Fenton  (1974)  for  w  >  0.92.  On  the  other 
hand  they  agree  with  the  numerical  results  of  Byatt-Smith  and  Longuet-Higgins 
(1976)  for  u  <  0.96.  Byatt-Smith  and  Longuet-Higgins  (1976)  also  used  an 
integro-differential  formulation*  However  they  were  not  able  to  compute  waves 
for  <d  >  0.96  because  too  many  mesh  points  were  required  to  describe 
accurately  the  flow  in  the  neighborhood  of  the  crest.  In  the  present  work 
this  difficulty  is  avoided  by  concentrating  the  mesh  points  near  the  crest  by 
an  appropriate  change  of  variable.  This  enables  us  to  compute  accurate 
solutions  up  to  u  -  0.99.  An  extrapolation  of  these  results  shows  that 
a  ♦  0.833  as  u  ♦  1.  This  constitutes  an  important  check  on  the  consistency 
of  our  two  numerical  schemas. 

The  problem  is  formulated  in  Section  2  and  the  highest  wave  is  calculated 
in  Section  3.  In  Section  4  we  compute  solitary  waves  of  arbitrary  amplitude 
via  an  integro-differential  equation  formulation.  The  results  are  discussed 
in  Section  5. 

2.  Formulation 

He  consider  a  two  dimensional  solitary  wave  in  an  inviscid  incompressible 
and  irrotational  fluid  bounded  below  by  an  horizontal  bottom,  we  take  a  frame 
of  reference  with  the  x-axis  parallel  to  the  bottom  and  moving  with  the  phase 
velocity  c  of  the  wave.  The  level  y  ■  0  is  chosen  as  the  undisturbed 
level  of  the  free  surface  and  gravity  is  assumed  to  act  in  the  negative  y 
direction. 

He  introduce  the  potential  function  4(x,y)  and  the  stream  function 
4(x,y).  without  loss  of  generality  we  choose  4*0  at  the  crest  and  4*0 
on  the  free  surface.  We  denote  by  Q  the  value  of  4  on -the  bottom.  Then 


the  undisturbed  depth  H  is  given  by 


H  ■  *  .  (2.1 ) 

c 

We  introduce  dimensionless  variables  by  taking  H  as  the  unit  length 
and  c  as  the  unit  velocity.  We  choose  the  complex  potential 

f  -  ♦  +  !♦  (2.2) 


as  the  independent  variable. 

We  shall  seek  the  complex  velocity 

C  -  u  -  iv  (2.3) 

as  an  analytic  function  of  f  in  the  strip  -1  <  ♦  <  0.  At  infinity  we 
require  the  velocity  to  be  c  in  the  x-direction  so  that  the  dimensionless 
velocity  is  unity  in  the  x-direction.  Therefore  C  must  tend  to  one  at 
infinity. 

On  the  free  surface,  the  Bernoulli  equation  yields 


Iu2(f)  +  v2(*)J 


♦ 


f#  v(») 

In  2  2 

W  u  (s)+v<e(s) 


F_ 

2 


-  a  on 


(2.4) 


Here  a  is  the  elevation  of  the  crest  and  F  is  the  Froude  number  defined  by 

F  -  —  .  (2.5) 

/gH 

The  functions  u(f)  and  v(+)  in  (2.4)  denote  respectively  u(+/°_)  and 
v(*,0_). 

On  the  bottom  the  kinematic  boundary  condition  yields 

v  -  0  on  ^  -  -1  .  (2.6) 

This  completes  the  formulation  of  the  problem  of  determining  the  analytic 
function  (.  This  function  must  tend  to  one  at  infinity  and  satisfy  (2.4)  on 
♦  *  1  and  (2.6)  on  *  ■  -1. 


Finally  let  us  mention  that  the  asymptotic  behavior  of  u( $)  -  iv(  $)  as 


♦  ♦  ±*  is  described  by  Stokes'  result 


u($)  -  iv($)  ~  Ae 


as  $  ♦  ±« 


(2.7) 


Here  A  is  a  complex  constant  to  be  found  as  part  of  the  solution  and  X  is 


the  smallest  root  of 


tan  wX 


(2.8) 


3.  The  highest  solitary  wave 

In  this  section  we  present  a  numerical  scheme  based  on  collocation  to 
compute  the  highest  solitary  wave.  This  wave  is  characterized  by  a  stagnation 
point  at  the  crest  where  the  surface  makes  a  120°  angle  with  itself.  (See 
Figure  1.) 

Following  Lenau  (1966)  we  introduce  the  new  variable  t  by  the  relation 


£  .  2  iog  i  . 


(3.1) 


This  transformation  maps  the  flow  domain  onto  the  domain  {|t|  <  1,  Im  t  >  0} 
in  the  complex  t  plane  (see  Figure  2).  He  use  the  notation  t  «  reio  so 
that  the  free  surface  is  described  by  r  -  1 ,  0  <  o  <  v. 

Lenau  (1966)  derived  the  following  expansion  for  the  complex  velocity  C 


5  -  «0(t) 


(3.2) 


where 


Q(t)  -  A(1-t2)2X  +  l  a  .tZn  . 


(3.3) 


Here  X  is  the  smallest  root  of  (2.8).  The  coefficients  A  and  a£  (i  * 
1,2,3,...)  in  (3.3)  have  to  be  found  to  satisfy  the  boundary  condition  (2.4), 
He  solve  the  problem  approximately  by  truncating  the  infinite  sum  in 
(3.3)  after  H  terms.  Differentiating  (2.4)  with  respect  to  o  and  using 


(3.1)  yields 


-5- 


F2[u(0)u  (o)  +  v(o)v  ( a)  ] 
0  0 


(3.4) 


*  ^2.  .  ~2  sin  0 
u  (0)+v  (0) 


■  0  . 


Here  u(0)  *  u[$(o)]  and  v(0)  -  v[$(o)]  are  the  components  of  the  velocity 
on  the  free  surface. 

«V  IV 

The  functions  u  and  v  and  their  derivatives  u  _  and  v„  are 

o  o 

io 

obtained  in  terms  of  A,  X  and  a^(i  “  1,...,N)  by  substituting  t  -  e  in 
(3.2).  We  find  the  N  +  3  unknowns  X,  A,  F,  a^(i“1 , . . . ,N)  by  satisfying 


(3.4)  at  the  N  +  2  mesh  points 


°I  "  2 (N+2 )  (I  "  2)  1-11 


■  N+2  . 


(3.5) 


Thus  we  obtain  a  system  of  N+2  nonlinear  algebraic  equations.  The  last 
equation  is  obtained  by  imposing  (2.8). 

This  system  of  M  +  3  equations  for  the  N  +  3  unknowns  X,  a,  F, 
a^(i««1 ,  • . . ,N)  was  solved  by  Newton's  iterations.  For  most  calculations  the 
values  a1  »  -  aA  -  0  (i*2,...,N),  A  -  -0.32,  X  -  0.32  and  F  -  1.3  were 
used  as  the  initial  guess.  The  method  converges  rapidly  and  a  residual  error 
of  10-10  was  obtained  after  4  or  5  iterations. 

Numerical  values  of  f  for  various  values  of  N  are  shown  in  Table  1 . 


These  results  indicate  that  the  value  F  *  1.29091  is  correct  to  5  decimal 

places.  The  profile  of  the  wave  is  shown  in  Figure  1 . 

Relation  (2.4)  shows  that  the  amplitude  a  of  the  highest  wave  is 

max 

given  by 


a  ■  -  0.83322  . 

max  2 


(3.6) 


This  value  is  about  0.006  higher  than  the  values  obtained  by  Yamada  (1957), 
Lenau  (1966),  Yamada,  et  al  (1968)  and  Longuet- Higgins  and  Fenton  (1974).  On 
the  other  hand  it  agrees  with  the  values  mentioned  by  Witting  (1975)  and 


Schwarts  and  Fenton  (1982) 


Our  numerical  Method  differ*  from  that  of  Lenau  (1966)  because  we  satisfy 


(3*5)  at  the  mesh  points  (3.6)  instead  of  solving  for  the  Fourier 
coefficients.  It  is  also  mo re  accurate  because  we  retain  up  to  100  terms  in 
(3.3)  whereas  Lenau  retained  only  9  terms. 

As  a  further  check  on  our  calculations  we  repeated  Yamada's  (1957) 
calculations.  Yamada  (1957)  presented  the  value  0.827  i  0.008  obtained  with 
11  mesh  points.  With  11  mesh  points  we  also  obtained  0.827.  However,  we 
obtain  0.832  with  30  mesh  points  and  0.833  with  100  mesh  points.  Thus 
Yamada's  (1957)  scheme  yields  the  correct  answer  when  a  sufficiently  large 
number  of  mesh  points  is  used. 


4.  numerical  solution  via  an  integro- differential  equation 

It  is  convenient  to  reformulate  the  problem  as  an  integro-differential 
equation  by  considering  u  -  iv  -  1.  This  function  tends  to  zero  at 
-infinity.  In  order  to  satisfy  the  boundary  condition  (2.6)  on  ♦  -  -1  we 
reflect  the  flow  inthe  boundary  ♦  -  -1 .  Thus  we  seek  u  -  iv  -  1  as  an 
analytic  function  of  f  in  the  strip  -2  <  ^  <  0. 

The  values  of  u  and  v  on  the  f ree  surface  -  0  and  its  image 

♦  »  -2  are  related  by  the  identities 

u(*,0)  -  u(*,-2) 
v($,0)  -  -v($,-2)  .  (4.2) 

In  order  to  find  a  relation  between  u(+,0)  and  v(  $,0)  we  apply 
Cauchy's  theorem  to  the  function  u  -  iv  -  1  in  the  strip  -2  <  ♦  <  0.  Using 
(4.1)  and  (4.2)  and  exploiting  the  bilateral  symmetry  of  the  wave  about 

♦  ■  0  we  obtain  after  some  algebra 


"«♦>  - 1  ■  p  £  *  ^Id* 


J 


■o 

.  4 


i 


i 


.1  r  i sJLHaljt  iMlkll  a.  (4.3) 

<a-4)  +  4 

+ 1  r  («+♦>*<«>  ±  2[u(«)-n  . 

«  /ft  2 

0  (s+4)2  +  4 

The  first  integral  in  (4.3)  is  of  Cauchy  principal  value  form.  We  shall 
measure  the  amplitude  of  the  wave  by  the  parameter  <a>  Using  the  symmetry  of 


the  wave  about  4*0/  we  rewrite  (1.2)  in  the  form 

w  »  1  -  F2[u(0)]2  . 

Using  (4.4)  and  (2.4)  evaluated  at  4*0  we  obtain 


(4.4) 


i 


i 


V 


+  .  (4.5) 

For  a  given  value  of  w,  (2.4),  (4.3)  and  (4.5)  define  a  system  of 
integro-differential  equatins  for  u(  4)  /  v(  4) ,  a  and  F. 

In  order  to  solve  these  equations  we  find  it  convenient  to  introduce  the 
new  variable  0  instead  of  4  by  the  relation 

4  -  0T,  Y  >  1  .  (4.6) 

Therefore  we  rewrite  (2.4),  (4.3)  and  (4.5)  in  terms  of  0,  u  ( 0)  -  u[4(0)] 
and  v*(0)  ■  v{ 4(0)1  • 

Next  we  introduce  the  N  mesh  points 

> ,M  (4.7) 

where  E  is  the  interval  of  discretization.  The  change  of  variable  (4.6)  is 

chosen  because  it  concentrates  the  mesh  points  near  the  crest  of  the  wave. 

For  very  steep  waves  the  value  of  y  was  taken  as  3. 

we  shall  satisfy  (2.4)  and  (4.3)  at  the  points  0  _ 

I+  2 

I  "  Thus  we  obtain,  after  discretization  2M  -  2  nonlinear 

•  * 

algebraic  equations  for  the  2M  +  2  unknowns  a,  F,  and  u  (0j),  v  (0  )  I  ■ 
1 , • • • ,M.  Relations  (4.4)  and  (4.5)  provide  two  more  equations.  An  extra 


0X  -  (1-1  )E  1*1, 


11 


(4.8) 


equation  is  obtained  by  imposing  the  symmetry  condition 

v  ( 0  )  *  0  . 

*  *  * 

The  last  equation  expresses  u  (0  )  in  term  of  u  (8  . )  and  u  ( 0  . )  by 

N  M— 1 

an  extrapolation  formula  based  on  the  asymptotic  formula  (2.7).  The 
discretization  of  (2.4)  and  (4.3)  follows  closely  the  work  of  Vanden-Broeck 
and  Schwartz  (1979). 

The  system  of  2M  +  2  equations  was  solved  by  Newton's  iterations. 

The  most  important  source  of  error  in  the  numerical  scheme  arises  from 

the  truncation  of  the  infinite  integrals  in  (4.3)  at 

a  "  ♦  -  t(M-1)E] Y  .  (4.8) 

max 

We  used  two  different  methods  to  approximate  the  infinite  integrals  in 

(4.3).  In  the  first  method  we  used  the  asymptotic  formula  (2.7)  to 

approximate  the  integrals  between  +  and  infinity.  This  approach  is 

max 

similar  to  the  method  used  by  Byatt- Smith  and  Longuet-Higgins  (1976).  In  the 

second  method  we  simply  neglected  the  contribution  of  the  integrals  between 

♦Mx  *nd  infinity.  In  this  second  method  we  also  replaced  the  equation  in 

which  an  extrapolation  based  on  (2.7)  is  used  by  a  Lagrange  extrapolation 

formula.  Thus  the  second  method  is  completely  independent  of  (2.7).  Both 

methods  were  found  to  give  accurate  results.  However  the  first  method  is  more 

efficient  because  accurate  results  can  be  obtained  with  $  relatively 

max 

small.  Most  of  the  results  presented  in  the  next  section  were  obtained  by 
using  the  first  method. 


5.  Discussion  of  the  results 

In  the  first  calculation  the  iterations  were  started  with  the  classical 

solution  of  the  Korteweg  and  de  Vries  equation.  For  o>  small  the  iterations 

converged  rapidly.  Once  a  solution  was  obtained  it  was  used  as  an  initial 

guess  for  a  larger  value  of  «  and  so  on. 

For  each  value  of  w  we  took  E  small  enough  and  $  large  enough 

max 

for  the  results  to  be  independent  of  E  and  $  .  This  was  achieved  in  the 

max 

following  way.  For  a  given  value  of  4  we  progressively  decreased  E  to 

max 

a  value  for  which  the  results  were  independent  of  E  to  the  degree  of 

accuracy  desired.  We  repeated  the  procedure  for  larger  and  larger  values  of 

$  up  to  a  value  for  which  the  results  were  also  independent  of  $ 
max  max 

This  procedure  is  illustrated  in  Table  2. 

In  Figure  3  we  present  the  numerical  values  of  the  Froude  number  F 
versus  u.  These  results  confirm  that  the  highest  solitary  wave  is  not  the 
fastest,  we  also  show  the  results  obtained  by  Longuet-Higgins  and  Fenton 
(1974)  and  by  Byatt-Saith  and  Longuet-Higgins  (1976).  Our  results  agree  with 
those  of  Longuet-Higgins  and  Fenton  for  u  <  0.92  and  with  those  of  Byatt- 
Smith  and  Longuet-Higgins  for  u  <  0.96. 

Byatt-Smith  and  Longuet-Higgins  (1976)  were  not  able  to  compute  waves  for 
u  >  0.96  because  their  numerical  procedure  uses  equal  increments  in  the 
velocity  potential.  This  is  not  well  suited  to  the  calculation  of  very  steep 
waves  because  large  curvature,  low  velocity  and  sparse  point  spacing  are 
characteristic  of  the  crest  region.  In  the  present  work  this  difficulty  has 
been  avoided  by  concentrating  the  mesh  points  near  the  crest  by  the  change  of 


variable  (4.6) 


Table  2s  Values  of  F  when 
u  -  0.98  and  y  •  3 


(b)  *  -  10 


F  N  F 


1.29141  100  1.29395 
1.29143  120  1.29152 
1.29145  150  1.29145 
1.29145  180  1.29145 


FIGURE  3 


Our  calculation*  show  that  tha  result*  of  Longuet-Higgins  and  Fenton 
(1974)  are  not  accurate  for  a  >  0.92.  This  seem*  to  confirm  Hitting'*  (1975) 
suggestion  that  their  method  is  defective  because  the  assumed  expansion  is 
incomplete. 
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Figure  1 «  Computed  free- surface  profile  for  the  highest  solitary  wave. 

The  vertical  scale  is  the  same  as  the  horizontal  scale. 

Figure  2.  Flow  configuration  in  the  complex  t-plane 

Figure  3.  The  Froude  number  F  as  a  function  of  at  as  given  by  the 

numerical  scheme  of  Section  4  (curve  a).  By a tt- Smith  and  Longuet 
Riggins  (1976)  (curve  b)  and  Longuet-Higgins  and  Fenton  (1974) 
(curve  c).  The  cross  corresponds  to  the  highest  wave  calculated 
in  Section  2. 
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